Determination of the characteristic curves of a nonlinear first order system from Fourier analysis

Based on Fourier analysis, we develop an expression for modeling and simulating nonlinear first order systems. This expression is associated to a nonlinear first order differential equation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=f(x)+g(x)x'$$\end{document}y=f(x)+g(x)x′, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x=x(t)$$\end{document}x=x(t) is the dynamical variable, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y=y(t)$$\end{document}y=y(t) is the driving force, and the f and g functions are the characteristic curves which are associated to dissipative and memory elements, respectively. The model is obtained from the sinusoidal response, specifically by calculating the Fourier analysis of y(t) for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x(t)=A_1\sin (\omega t)+A_0$$\end{document}x(t)=A1sin(ωt)+A0, where the model parameters are the Fourier coefficients of the response, and the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_0$$\end{document}A0, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_1$$\end{document}A1 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A_1'=A_1\omega$$\end{document}A1′=A1ω. The same expression is used for two kinds of time-domain simulations: to calculate other driving force \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{y}$$\end{document}y^ based on a dynamical variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{x}$$\end{document}x^; and, to calculate the dynamical variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{x}$$\end{document}x^ based on a driving force \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{y}$$\end{document}y^. In both cases, the dynamical variable must remain in the range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{x}\in [A_0-A_1,A_0+A_1]$$\end{document}x^∈[A0-A1,A0+A1]. By analyzing this expression, we found an equivalence between the Fourier coefficients and the polynomial regressions of the characteristic curves of f and g. This result allows us to obtain the system modeling and simulation based on the amplitude and phase Fourier spectrum obtained from the Fast Fourier Transform (FFT) of the sampled \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_n$$\end{document}yn version of y(t). It is shown that this technique has a low computational complexity, and it is expected to be suitable for real-time applications for system modeling, simulation and control, in particular when the explicit expressions of the characteristic curves are unknown. Fourier analysis is a fundamental tool in electronics, mathematics and physics, but to the best of the author’s knowledge, no work has found this clear evidence of the connection between the Fourier analysis and a first order differential equation. The aim of this work is to initiate a systematic study on this topic.


Supplementary M1. Relation between FFT and Fourier Series
The Fast Fourier Transform (FFT) is an algorithm that calculates the Discrete Fourier Transform (DFT) of a sequence. The DFT transforms a sequence of N complex numbers {x n } := {x 0 , x 1 , · · · , x N−1 } into another sequence of complex numbers {X k } := {X 0 , X 1 , · · · , X N−1 }, which is defined by 1 with k ∈ [0, N − 1]. The Inverse Discrete Fourier Transform (IDFT) is obtained by with n ∈ [0, N − 1]. Consider a real function x(t) with period T and t ∈ R, it can be expressed as a Fourier Series (FS) (a k cos(kωt) + b k sin(kωt)) .
If x(t) is sampled with period T s , its sampled sequence x n = x(t n = nT s ) and the sequence length N is chosen to verify x(t 0 ) = x(t 0 + NT s ), then by using the Whittaker-Nyquist-Shannon sampling theorem 2 , the original signal x(t) can be fully reconstructed by the expression x n =x(t n ) = x(nT s ) = a 0 + where ⌊⌋ is the floor function. Note that the maximum frequency is where T = NT s and f s = 1/T s is the sampling frequency. The last inequality can be obtained by separating the cases Equation (5) is consistent with the sampling Shannon condition f s > 2 f max for a fully reconstructed signal. In the following, we will find the relation between the Fourier coefficients of Eq. (4) with those of Eq. (2). For real sequences x n it is verified , where * denotes the complex conjugate, and ⌊ ⌋ is the floor function. This result follows directly from Eq. (1) by using the relation x n = x * n , which is valid for real sequences. By using Eq. (7), we can expand Eq. (2) as 2Re(X k e i 2π N kn ) + (−1) n (1 + (−1) N ) 2N By replacing X k = Re(X k ) + iIm(X k ) into Eq. (8) we obtain Finally, by comparing Eqs. (9) and (4) we identify Most of the numerical FFT methods calculate the double-sided DFT, which is consistent with the definition of Eq. (1). In order to calculate single-sided DFT, we use Eqs. (10) and (11).
A Matlab script to obtain the single-sided FS coefficients from the FFT function is presented here

Supplementary M2. Complex dynamical variable
The equations which were deduced in this work are based on a real dynamical variable x(t), which is forced to be equal to x(t) = A 1 sin ωt + A 0 for the system modeling. This appendix considers an extension to the complex domain, where the dynamical variable is a complex functionx(t) . For the system modeling, consider that it verifiesx(t) = re iωt +Ã 0 , with r > 0 andÃ 0 a complex number. As it is shown below, the mathematical expressions which arise in this complex formalism are simpler compared to the real formalism. However, the real formalism which is used throughout this work can be more directly applied to experimental or simulation data. Suppose thatỹ(t) is a complex-valued function of a real variable t with period T , we can expand it with the complex Fourier series where ω = 2π/T is the fundamental frequency and By replacing C k = Re(C k ) + iIm(C k ) into Eq. (12), we obtaiñ By usingx(t) = re iωt +Ã 0 andx ′ (t) = iωre iωt , Eq. (14) can be rewritten as whereÃ ′ 1 :=Ã 1 ω. In this complex version, the polynomial expansion analog to Eqs. (16) and (17) of the main text are in fact decoupled, and can be directly identified as In the following, we consider the discrete version of the complex formalism, which may be useful for a practical implementation based on a complex sampled data {x n ,ỹ n } for n ∈ [0, N − 1]. As a result of this analysis, we obtain the relation between the complex Fourier series and the FFT of the complex functionỹ n .

Supplementary M3. Computational complexity
On the one hand, the model is obtained from FFT, which has a computational complexity of O(N ln N) 1 , where N is the length of the sequence y n . On the other hand, the simulation by using Eq.

Supplementary M5. Demonstrative examples
This section presents some demonstrative examples of the formalism of this work. We first analyze a nonlinear inductor in series with a discontinuous resistance in order to show the complete procedure. Then we consider more simple and illustrative examples in order to discuss some specific details of the formalism. We will follow the steps that are presented in figure 1 of the main text.

4/12
Nonlinear inductor in series with a discontinuous resistance In this example, we consider a nonlinear RL system with a discontinuity, where the R and L components are defined by According to the section Nonlinear series RL system of the main text, the dynamical variable for this system is the current i(t), the driving force is the voltage v(t), and the system satisfies the first order differential equation The characteristic curves are f (i) = R(i) · i and g(i) = L(i). In order to obtain the system modeling, we use a dynamical variable i(t) = A 1 sin(ωt), with A 1 = 4 A and ω = 2π50 rad/seg. Hz, which corresponds with the 40th harmonic, is at least 100 times smaller than the first harmonic. This indicates that a model with 40 harmonics is probably enough to fully capture the dynamics. Then, the Fourier coefficients can be calculated by using where P k is the single-sided FFT of the signal and the sub-index identifies the number of harmonic, which in this example will be constrained to k ∈ [0, 40]. Another equivalent form to obtain the Fourier coefficients is from the real and imaginary parts of the FFT single-sided spectrum, by using Eqs. (10) and (11) from Supplementary M1, see also the Matlab script at the end of that Supplementary. In this example we obtain a total of 2 · 40 + 1 = 81 Fourier coefficients for the system modeling. It is worth mentioning that this simulations is noiseless, but in case of noise presence, use can be made of the wide range of techniques of spectral density estimation 3 in order to obtain a cleaner Fourier spectrum. The next step is Step ( figure 5 shows an oscillatory effect at the discontinuity that is typical of a Gibbs-Wilbraham phenomenon [4][5][6] . In order to quantitatively assess the error between the theoretical simulation and the model, we compute the root-mean-square error (RMSE). We obtain RMSE = 9.66 for the f function and RMSE = 1 × 10 −5 for the g function. The g curve has been correctly identified and there is no presence of the discontinuity, this is a consequence of that f and g depends on different Fourier coefficients, as shown in Eqs. At this stage, we have estimated the characteristic curves f and g. In order to test the model, we consider a simulation according to Step (4) in figure 1 of the main text. In a first test, we address the kind (i) of simulation by setting a dynamical variableî(t) = sin(2π50 t) + sin(2π100  figure 11. This figure shows more clearly the presence of the Gibbs phenomenon at the discontinuities. A second test is presented in the following in order to visualize the kind (ii) of simulation. We set a driving forcê v(t) = 300 sin(2π50 t) + 90 and calculate the corresponding dynamical variableî, with the initial condition i(0) = 2.5 A. We use Eq. (38) of the main text to obtain the corresponding dynamical variable. The results are shown in supplementary figure 12, where the solid line corresponds to the theoretical simulation of Eq. (23) and the circles are the predicted values of the model by using Eq. (38). The RMSE is 0.049, where the mayor discrepancy is around the discontinuity at i = 1A, this is attributed to the Gibbs phenomenon mentioned above. The dataset for the simulation is obtained with a simulation step T s = 10 µseg and a length N = 4000. The predicted values by using Eq. (38) is plotted every 100 µseg for ease of visualization, this means that there is 400 circles in the figure.

7/12
Diode An ideal diode is one of the simplest devices which can be defined by a nonlinear I-V curve. In this section, we apply the RC parallel system to model a diode (see the section Nonlinear parallel RC system of the main text). The nonlinear R component represents the static nonlinear I-V curve of the diode, and the C component can be associated to a loss capacitor. Consider the dynamical variable v(t) = A 1 sin(ωt), with A 1 = 0.8 V and ω = 2π50 rad/seg. The temporal evolution of this dynamical variable and the driving force, which corresponds to the diode current, are shown in supplementary figure 13. The FFT of the driving force is shown in supplementary figure 14. It can be seen that the amplitude is negligible at 1000 Hz, which corresponds to the 20th harmonic, however we will choose up to the 40th harmonics consistent with the previous example. Now, we use Eq. (36) of the main text in order to calculate the characteristic curves. In this example, the g function, which corresponds to the capacitor has neglected values. For this reason, we decide to calculate the I-V curve directly, this is shown in supplementary figure 15. We obtain RMSE(i(t) −î(t)) = 2 · 10 −5 . We test the model with a different dynamical variablev(t) = 0.3 sin(2π50 t) + 0.3 sin(2π100 t) + 0.3 sin(2π150 t), and as a result, we obtain an error RMSE(i(t) −î(t)) = 6 · 10 −6 . For a chirp dynamical variable, with 0.8 V in amplitude and an increasing frequency from 0 to 500 Hz in a simulation time of 50 mseg, we obtain RMSE(i(t) −î(t)) = 2 · 10 −5 .

8/12
Diode in parallel with a capacitor Consider a circuit with a diode in parallel with a 100 µF capacitor. To model the system, we apply a dynamical variable v(t) = A 1 sin(ωt), with A 1 = 0.8 V and ω = 2π50 rad/seg. This dynamical variable and the corresponding driving force are shown in supplementary figure 16. The FFT of the driving force is shown in supplementary figure 17. Then, we use Eq. (36) of the main text to verify the prediction of the I-V characteristic curve, this is shown in supplementary figure 18. We obtain RMSE(i(t) −î(t)) = 1.1 · 10 −5 . For a different dynamical variable,v(t) = 0.3 sin(2π50 t) + 0.3 sin(2π100 t) + 0.3 sin(2π150 t) we use the model to compute the driving force, obtaining RMSE(i(t) −î(t)) = 2.3 · 10 −5 , where the correspondingly I-V curve is shown in supplementary figure 19. For a chirp dynamical variable, with 0.8 V in amplitude and an increasing frequency from 0 to 500 Hz, in a simulation time of 50 m seg, we obtain RMSE(i(t) −î(t)) = 2.5 · 10 −4 , where the correspondingly I-V curve is shown in supplementary figure 20. This example shows that the estimation of the characteristic curves from the sinusoidal response allows us to obtain a model which is able to predict the response from a wide variety of new dynamical variables, in particular a chirp response. Notice that this has been possible because the system follows a first order differential equation such as that of Eq. (44) of the main text. However, many real systems are frequency dependent, this means that the resistive and inductive elements are explicitly frequency dependent. In these cases, the formalism of this work is expected to be valid only on a frequency range around the frequency that was used for the system modeling. See the Discussion section in the main text for more details.

Nonlinear inductor
In the main text, we have shown that the Fourier analysis of the driving force is equivalent to the polynomial regressions of the characteristic curves, and the polynomial regressions are known to present the Runge phenomenon. The aim of this example is to show the presence of the Runge phenomenon. This section studies a nonlinear inductor in series with a resistance of 0.1 Ω. We use the constitutive equation for the inductor where L(i) = 0.02 1 − tanh 2 (i/5) .
According to Step (1) in figure 1 of the main text, we use a sinusoidal dynamical variable for the system modeling, in particular, we set i(t) = A 1 sin(ωt), with A 1 = 25 A and ω = 2π50 rad/seg. The dynamical variable and the driving force are shown in supplementary figure 21. We compute the FFT of the driving force according to Step (2), this is shown in supplementary figure 22. It can be seen that the amplitude is almost negligible at 1500 Hz, which corresponds to the 30th harmonic, then, we consider a system modeling up to the 40th harmonic. We test the model with the same dynamical variable by using Eq. (36) of the main text, this is shown in supplementary figure 23. The obtained error is RMSE(v(t) −v(t)) = 0.37. Supplementary figure 24 shows the characteristic curve L(i) from Eq. (26) compared to the g function calculated by Eq. (36) of the main text. An oscillatory Runge phenomenon 7 , which is typical of a polynomial curve fitting, can be observed at the extreme values. This oscillation can be improved by reducing the order k max , but at the expense of a worst fit. Now, once the characteristic curves have been obtained, the model is ready to be tested. In this example, we use the kind (i) of simulation according to figure 1 of the main text. In the following we apply dynamical variables restricted to |i(t)| < 20 A in order to avoid the range of values where the Runge oscillation is present. Firstly, we consider the sum of two sinusoidal functionsî(t) = 10 sin(ωt) + 10 sin (